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We study in detail the hydrodynamic theories describing the transition to collective motion in 
polar active matter, exemplified by the Vicsek and active Ising models. Using a simple phenomeno¬ 
logical theory, we show the existence of an infinity of propagative solutions, describing both phase 
and microphase separation, that we fully characterize. We also show that the same results hold 
specifically in the hydrodynamic equations derived in the literature for the active Ising model and 
for a simplified version of the Vicsek model. We then study numerically the linear stability of these 
solutions. We show that stable ones constitute only a small fraction of them, which however in¬ 
cludes all existing types. We further argue that in practice, a coarsening mechanism leads towards 
phase-separated solutions. Finally, we construct the phase diagrams of the hydrodynamic equations 
proposed to qualitatively describe the Vicsek and active Ising models and connect our results to the 
phenomenology of the corresponding microscopic models. 
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I. INTRODUCTION 

Collective motion is the ability of large groups of motile 
agents to move coherently on scales much larger than 
their individual sizes. It is encountered at all scales in 
nature, from macroscopic animal groups, such as bird 
flocks, fish schools, or sheep herds, down to the cellu¬ 
lar scale, where the collective migration of cells |T] or 
bacteria [3] is commonly observed. At the subcellular 
level, in vitro motility assays of actin filaments |5] or mi¬ 
crotubules |T] have shown the spontaneous emergence of 
large vortices. Collective motion is also observed in en¬ 
semble of man-made motile particles such as shaken polar 
grains |5], colloidal rollers |B], self-propelled droplets 0 
or assemblies of polymers and molecular motors 1311E]. 
Despite the differences in their propulsion and interac¬ 
tion mechanisms, these seemingly very different systems 
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FIG. 1. Top: Micro-phase separation in the Vicsek model, rj = 0.4, vq = 0.5, po = 0.83, 1.05, 1.93. Bottom: Phase separation 
in the Active Ising model. D = 1, e = 0.9, /3 = 1.9, po = 1.5, 2.35, 4.7. System sizes 800 x 100. High-density bands propagate 
as indicated by the red arrows on the left snashots. 


share common macroscopic behaviors that can be cap¬ 
tured by minimal physical models. Of particular inter¬ 
est is the emergence of directed collective motion, which 
was first addressed in this context in a seminal work by 
Vicsek and co-workers [H]. The Vicsek model consists 
in point particles moving at constant speed and aligning 
imperfectly with the direction of motion of their neigh¬ 
bors. When the error on the alignment interaction is de¬ 
creased, or the density of particles increased, a genuine 
phase transition from a disordered to a symmetry-broken 
state is observed. This flocking transition gives rise to an 
emergent ordered phase, with true long-range polar or¬ 
der even in 2D, where all the particles propel on average 
along the same direction. Toner and Tu showed ana¬ 
lytically, using a phenomenological fluctuating hydrody¬ 
namic description, how this ordered state, which would 
be forbidden by the Mermin-Wagner theorem at equi¬ 
librium [in], is stabilized by self-propulsion m- The 
transition to collective motion in the Vicsek model has a 
richer phenomenology than originally thought. As first 
pointed out numerically in m, at the onset of collec¬ 
tive motion, translational symmetry is broken as well. 
In periodic simulation boxes, high-density ordered bands 
of particles move coherently through a low-density disor¬ 
dered background. The transition between these bands 
and the homogeneous disordered profile is discontinuous, 
with metastability and hysteresis loops. These spatial 
patterns and the first-order nature of the transition can 
be encompassed in a wider framework, which describes 
the emergence to collective motion as a liquid-gas phase 
separation nsun]. The travelling bands result from the 
phase coexistence between a disordered gas and an or¬ 
dered polarized liquid. This framework captures many of 
the characteristics of the transition, from the scaling of 
the order parameter to the shape of the phase diagram. 
This phase-separation picture is robust to the very de¬ 
tails of the propulsion and interaction mechanisms. More 


specifically, it has also been quantitatively demonstrated 
in the active Ising model m in which particles can dif¬ 
fuse in a 2d space but self-propel, and align, only along 
one axis. However the specifics of the emergent spatial 
patterns and the type of phase separation depend on the 
symmetry of the orientational degrees of freedom. While 
the active Ising model model shows a bulk phase separa¬ 
tion, the Vicsek model is akin to an active XY model and 
is associated with a microphase separation where the co¬ 
herently moving polar patterns self-organize into smectic 
structures HD (see Fig.0. 


In this paper, building on the two prototypical models 
that are the Vicsek model and the active Ising model, we 
provide a comprehensive description of the emergent pat¬ 
terns found at the onset of the flocking transition from a 
hydrodynamic perspective. We first recall the definitions 
and phenomenologies of these two models in Sec. |TT| In 
Sec. |III| we provide a simplified hydrodynamic descrip¬ 
tion of the flocking models. In line with US! US], we show 
that these models support non-linear propagative solu¬ 
tions whose shape is described using a mapping onto the 
trajectories of point-like particles in one-dimensional po¬ 
tentials. Finding such solutions thus reduces to a classical 
mechanics problem with one degree of freedom. For given 
values of all the hydrodynamic coefficients, and hence of 
all underlying microscopic parameters, we And an infinity 
of solutions, describing both phase and microphase sepa¬ 
rations, that we fully characterize. We then show that the 
same results hold specifically for the hydrodynamic equa¬ 
tions explicitly derived for the active Ising model m and 
for a simplified version of the Vicsek model m- Next, 
we investigate the linear stability of these solutions as 
solutions of the hydrodynamic equations in Sec. |VI| and 
their coarsening dynamics in Sec. VII Finally, we pro¬ 
vide full phase diagrams constructed from the hydrody¬ 
namic model in Sec. |VIII[ We close by discussing the 
similarities and differences with the phenomenology of 
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the agent-based models and conjecture on the role of the 
hydrodynamic noise in the selection of the band patterns. 


II. PHENOMENOLOGY OF MICROSCOPIC 
MODELS 


Let us first briefly recall the phenomenology of the Vic- 
sek and active Ising models. They are both based on the 
same two ingredients: Self-propulsion and a local align¬ 
ment rule. The major differences between the two models 
are thus the symmetries of the alignment interaction and 
of the direction of motion. 



FIG. 2. Phase diagrams of the microscopic models. The 
red and blue lines delimit the domain of existence of (micro) 
phase-separated profiles. The black line and squares indicate 
the position of the snapshots shown in Fig. vg = 0.5 for the 
Vicsek model, D = 1 and e = 0.9 for the active Ising model. 


A. Vicsek model 

In the Vicsek model [5] , JV point-like particles, labeled 
by an index i, move at constant speed Vg on a rectan¬ 
gular plane with periodic boundary conditions. At each 
discrete time step At = 1, the headings di of all particles 
are updated in parallel according to 

di(t+ 1) = {0j{t))j^_\f. +r]^l (1) 

where A/) is the disk of unit radius around particle i, a 
random angle drawn uniformly in [—tt, tt], and r] sets the 
level of noise, playing a role akin to that of a temperature 
in a ferromagnetic XY model. Then, particles hop along 
their new headings: r.i(t-l-l) = ri(t)where 
is the unit vector pointing in direction given by Oi(t + 1). 


B. Active Ising model 

In the active Ising model m. particles carry a spin 
±1 and move on a 2D lattice with periodic boundary 
conditions. Their dynamics depend on the sign of their 
spin: A particle with spin s jumps to the site on its 
right at rate D{1 + se) and to the site on its left at rate 
D{1 — se), where 0 < £ < 1 measures the bias on the 
diffusion. On average, -bl particles thus self-propel to 
the right and —1 particles to the left at a mean velocity 
Vg = 2De. Both types of particles diffuse symetrically at 
rate D in the vertical direction. 

The alignment interaction is purely local. On a site i, 
a particle flips its spin s at rate 

Wi{s-)--s)=exp(^-^'^'^ ( 2 ) 

where T is a temperature, and mi and pi are the magne¬ 
tization and number of particles on site i. (An arbitrary 
number of particles is allowed on each site since there is 
no excluded volume interaction.) 


C. A liquid-gas phase transition 

The phase diagrams in the temperature/noise-density 
ensemble are shown for both models in Fig. high¬ 
lighting their similarity. At high temperature/noise or 
low density both systems are in a homogeneous disor¬ 
dered gas state. At low temperature/noise and high den¬ 
sity they are homogeneous and ordered; in these liquid 
phases, all particles move in average in the same direc¬ 
tion. In the central region of the phase diagram, inhomo¬ 
geneous profiles are observed, with liquid domains mov¬ 
ing in a disordered gaseous background. 

The phase transitions of both models have all the fea¬ 
tures of a liquid-gas transition, exhibiting metastability 
and hysteresis close to the transition lines [IMl]- The 
main difference between the two models lies in the co¬ 
existence region: In the active Ising model, the particles 
phase separate in a gaseous background and an ordered 
liquid band, both of macroscopic sizes m- The coexist¬ 
ing densities depend only on temperature and bias, but 
not on the average density; in the coexistence region, in¬ 
creasing the density at fixed T, e thus results in larger and 
larger liquid domains whose density remains constant, as 
shown in Fig. [l] Conversely, in the Vicsek model, the 
system forms arrays of ordered bands arranged periodi¬ 
cally in space which have a finite width along their direc¬ 
tion of motion: a micro-phase separation occurs |14j . As 
shown in Fig. increasing the density at constant noise, 
the number of bands increases but their shape does not 
change El- 

Three types of propagating patterns can thus be ob¬ 
served at phase coexistence, all shown in Fig.[^ (i) local¬ 
ized compact excitations, (ii) Smectic microphases, and 
(iii) Phase-separated polar liquid domains. In the vicin¬ 
ity of the left coexistence line, collective motion emerges 
in the form of localized compact excitations in both mod¬ 
els [22]. At higher density, phase-separated domains are 
found in the active Ising model and periodic “smectic” 
bands in the Vicsek model. Understanding the emer¬ 
gence of these three types of solutions will be the focus 
of the rest of the paper. 
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III. HYDRODYNAMIC EQUATIONS 


corresponding to the Ising symmetry—are given by 


A lot of attention has been given in the literature to 
hydrodynamic equations of flocking models. Two differ¬ 
ent approaches have been followed, starting from phe¬ 
nomenological equations [niiiniin] or deriving explicitly 
coarse-grained equations from a microscopic model [HI 
HSl HHHSO]. All these equations describe the dynamics 
of a conserved density field p(f, t) coupled to a non- 
conserved magnetization field, the latter being a vector 
rn(r, t) for continuous rotational symmetries, as in the 
Vicsek model, or a scalar m(r,t), for discrete symme¬ 
tries, as in the active Ising model. 

We first introduce in Sec. |III A| two sets of hydrody¬ 
namic equations derived by coarse-graining microscopic 
models which will be discussed in this paper. Then, we 
turn in Sec. |IIIB| to a simpler set of phenomenological 
hydrodynamic equations on which we will establish our 
general results in Sec. [Tvl 


A. Coarse-grained hydrodynamic descriptions 


We first consider the equations proposed by Bertin 
et al. to describe a simplified version of the Vicsek 
model in which one solely considers binary colli¬ 
sions between the particles. One can then use, assum¬ 
ing molecular chaos, a Boltzmann equation formalism to 
arrive at the following hydrodynamic equations for the 
density field and a vectorial magnetization field [ 23 ] 

^ = -voV ■ rh (3) 

-I- 7 (to • V)m = — ^Wp+ ^V(|toP) 

— k{\ 7 • rh)m + {p — (^\rn^\)m (4) 

The mass-conservation equation (|^ simply describes the 
advection of the density by the magnetization field. 
Equation Q can be seen as a Navier-Stokes equation 
complemented by a Ginzburg-Landau term (/i—C|mp)TO, 
stemming from some underlying alignment mechanism, 
and leading to the emergence of a spontaneous magne¬ 
tization. Because particles are self-propelled in a given 
frame of reference, these equations break Galilean invari¬ 
ance so that one can have 7 yf 1 and k ^ 0 unlike e.g. in 
the Navier-Stokes equation. 

In Eqs. ^ and Q, to which we refer to as “Vicsek 
hydrodynamic equations” hereafter, all the coefficients 
7 , V, K, p and depend on the local density; see [T5| for 
their exact expression. 

The second set of equations, which we refer to as 
“Ising hydrodynamic equations” in the following, has 
been derived to describe the large-scale phenomenology 
of the active Ising model m- In this case, the dynam¬ 
ics of the density field and the scalar magnetization— 


— = DAp - vodxm (5) 

^ = DA'm-VQdxP + 2(^P ( 6 ) 

where /3 = 1 /T, a and r are positive coefficients depend¬ 
ing on /3 only and vq = 2De. The advection term VQ^Irh 
of Eq. ^ is here replaced by a partial derivative vodxm 
because, in the active Ising model, the density is advected 
by the magnetization only in the cc-direction. 


B. Phenomenological hydrodynamic equations 
with constant coefficients 


Goarse-grained hydrodynamic equations derived from 
microscopic models have the advantage of expressing the 
macroscopic transport coefficients in terms of microscopic 
quantities (noise, self-propulsion speed, etc). However, 
these possibly complicated relations may not be relevant 
to understand the qualitative behavior of the models. 
Thus, before discussing the Vicsek and Ising hydrody¬ 
namic equations in Sec. |V] we first study in detail, in 
Sec. jlVj a simpler model 

dtp = -vqV ■ rh (7) 

dtfh + ^{m ■ V)m = DV^m — AVp -I- a 2 m — a 4 |mpm 

( 8 ) 


where the transport coefficients vq, D, A and 04 are 
constant. In the following, we refer to these equations as 
the “phenomenological hydrodynamic equations”. This 
simplified model is very similar to that first introduced by 
Toner and Tu from symmetry considerations m- How¬ 
ever, unlike the original Toner and Tu model, we keep 
an explicit density dependence in 02 : 02 (p) = p — ipg, 
which is essential to account for inhomogeneous pro¬ 
files [inifnifTT] . 

The stability criteria of the homogeneous solutions 
(p(r,f) = po, m{r,t) = mg) of Eqs. ( [7p l are readily 
computed: 

• For Po < <Pg {^ 2 {po) < 0) only the disordered solu¬ 
tion (po, I Too I = 0 ) exists and is stable 


• For Po > Pg (02 (po) > 0) the disordered solution 
becomes unstable and ordered solutions (po, |mo| = 
a/(Po - Pg)/a4) appear. 

• The ordered solutions are linearly stable only when 

P0>Pe = Pg+ 4aJo+2X 

Thus, in the range po € [pg,pi], homogeneous solutions 
are linearly unstable. In the language of the liquid-gas 
transition, pg and pi are the gas and liquid spinodals, 
between which the homogeneous phases are linearly un¬ 
stable and spinodal decomposition takes place. In the 
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next section we address the existence of heterogenous or¬ 
dered excitations propagating through stable disordered 
(gaseous) backgrounds. This analysis will make it pos¬ 
sible both to identify all the possible flocking patterns 
and to further understand the first-order nature of the 
flocking transition. 



IV. PROPAGATIVE SOLUTIONS 

Let us now establish and classify all the inhomogeneous 
propagating solutions of Eqs. (7|8 1 . In order to do so, we 
first recast this problem into a dynamical system frame¬ 
work in section IIV Al We then show in section IIVBI that 
three types of propagating solutions exist with different 
celerities c and densities of the gaseous background pg. 


Sections IV C IV D and IV E are dedicated to a detailed 


study of how these solutions depend on c and Pg. Sec- 
tion |IV F1 shows how, once the average density is fixed, we 
are left with a one-parameter family of solutions. Last, 
section [TV G| is devoted to cases where the inhomogeneous 
profiles can be studied analytically. 


A. Newton Mapping 

Following m, we look for inhomogeneous polar exci¬ 
tations invariant along, say the y direction, and which 
propagate, and/or relax solely along the x direction. We 
thus assume, my = 0 and reduce Eqs. to: 


dtp = -vod^m ( 9 ) 

dtm + ^mdxm = Dd^m — XdxP + {p — <Pg)m — 04771 ^ 

( 10 ) 

where we wrote m = to ease the notation. We look for 
solutions propagating steadily at a speed c. Introducing 
the position z = x — ct in the frame moving at c: p{x, t) = 
p{z), m{x,t) = m{z), we obtain 

cp — vom = 0 (11) 

Dm -I- (c — ^m)m — Xp + {p — <Pg)m — a/pm^ = 0 (12) 

where the dots denote derivation with respect to z. Solv¬ 
ing Eq. (11 1 gives p{z) = Pg + ^m{z). If p{z) is localized 
in space, pg has a simple meaning. Since p{z) = pg when 
m{z) = 0 , the integration constant pg is the density in the 
gaseous phase surrounding the localized polar excitation. 
We can then insert the expression of p in Eq. (12 1 and 


obtain the second-order ordinary differential equation 

Dm+ (c-— ^m)rh— {ipg — pg)m-\ — -m?' — apm^ = 0 

(13) 

Following USEE], we now provide a mechanical inter¬ 
pretation of Equation (131 through the well-known New- 


FIG. 3. Left: Density and magnetization profiles of a prop¬ 
agative solution of the hydrodynamic Eqs. (|7]|^. Center: 
Magnetization profile in the comoving frame z = x — ct or, 
equivalently, trajectory m[z) of a point particle in the spuri¬ 
ous time a. Right: Phase portrait corresponding to the tra¬ 
jectory m{z). Parameters: D = vo = X = ^ = = (pg = 1. 




FIG. 4. The green potential can give rise to physical (positive, 
non-exploding) solutions while the red ones are ruled out by 
our conditions (Gl) (left) and (G2) (center). 


ton mapping. Rewriting Eq. (131 as: 


dH 


dm 


m? 

Vo : 
3c”^ 


f{m) = c — 


(14) 

(15) 

(16) 


it is clear that this equation corresponds to the mechan¬ 
ical equation of motion of a point particle. The position 
of the particle is m, z is the time variable, D is the mass 
of the particle, H{m) is an energy potential and f(m) 
is a position-dependent friction. The trajectory m{z) of 
this Active particle exactly corresponds to the shape of 
the propagative excitations of our hydrodynamic model 
in the frame moving at a speed c (see Fig. [^. 

We shall stress that for a given hydrodynamic model, 
Eq. (141 is parametrized by the two unknown parameters 


c and Pg which a priori can take any value. Each pair (c, 
Pg) gives different potential H and friction /, and hence 
different trajectories m{z). We now turn to the study 
of these trajectories and of the corresponding admissible 
values for the celerity c and the gas density Pg. 


B. Three possible propagating patterns 

The original problem of finding all the inhomogeneous 
propagative solutions m{x,t), p{x,t) of the hydrody¬ 
namic equations is now reduced to finding all the pairs (c. 
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FIG. 5. Example of the three types of trajectories. From left to right: magnetization and density profiles, phase portrait, and 
potential H. Top row: Periodic trajectory, pg = 0.835, c = 1.14. Center row: Homoclinic trajectory, pg = 0.83412, c = 1.14. 
Bottom row: Heteroclinic trajectory, pg = 0.83333, c = 1.1547. Phase portrait: Crosses indicate saddle points at m = 0 and 
m = m+. Squares indicate stable fixed points at m = m_. Potentials: The blue dashed line indicate where the friction changes 
sign. The red portion of the potential is the one visited by the trajectory. Parameters: D = vo = A = ^ = 04 = (/ 3 g = 1. 


Pg) for which the corresponding trajectories m{z) exist. 
Mass conservation, Eq. (§, imposes the boundary con¬ 
dition m{z = — 00 ) = m{z = -bcx)) so that we are looking 
for solutions of ( |14[ ) which are closed in the (m, m) plane. 
An example of propagative solutions m{x,t), p{x,t) to¬ 
gether with the corresponding trajectory m{z) and its 
phase portraits is shown in Fig. 

To put a first constraint on pg and c, let us rule out 
the potentials which cannot give such physical solutions. 
The extrema of H, solutions of H'{m) = 0, are located 
at m = 0 and m = m± with 


We can already discard the case where H'{jn) has two 
complex roots since H then has a single maximum at 
m = 0, and all trajectories wander to m = ±00 (see 
Fig. left). This leads to a first condition on c, pg-. 

- Pg)c^ < a^vl (Cl) 

Without loss of generality we can assume that c > 0 
and look only for solutions with m > 0. This rules out the 


(c, Pg) values for which m_ < 0 and m+ > 0 which give 
oscillations between negative and positive values of m 
(see Fig.|^ central panel). At the hydrodynamic equation 
level, such solutions would indeed correspond to different 
parts of the profiles moving in opposite direction. The 
corresponding condition 


Pg < ‘Pg 


(C2) 


imposes 0 < to_ < m+. The potential H then has two 
maxima, at m = 0 and m = rn+, and one minimum, at 
m = rri-. The typical shape of potential which gives ad¬ 
missible solutions is shown in Fig. along with examples 
of potentials ruled out by conditions (jcTl and (C21. 


>From the admissible shape of the potential H, we 
can now list all possible trajectories m(z) and the corre¬ 
sponding fields m{x,t), p{x,t)\ 


• Limit cycles, whose corresponding magnetisation 
profiles are periodic bands, as shown in the first 
row of Fig. 

• Homoclinic orbits, that start infinitely close to a 
maximum of H, hence spending an arbitrary large 
time there, before crossing twice the potential well 
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in a finite time to finally return to the same maxi¬ 
mum oi H at z = oo. These trajectories correspond 
to isolated solitonic band profiles, as shown in the 
second row of Fig. 

• Heteroclinic orbits that spend an arbitrary large 
time close to a first maximum of H, cross the po¬ 
tential well in a finite time, spend an arbitrary large 
time close to the second maximum of H, before re¬ 
turning to the first maximum. These trajectories 
correspond to phase separated profiles. The arbi¬ 
trary waiting times at the two maxima of H then 
reflect the arbitrary sizes of two phase-separated 
domains (see the third row of Fig. [^. 

A third condition on Pg,c arises from the non-linear 
friction term. Following the classical mechanics analogy, 
we define an energy function 


E = -Drri^ 
2 


H 


(18) 


Multiplying the equation of motion (131 by rh, we get 


dE 

dz 


—f{m)rh^ 


(19) 


Energy is injected when /(to) < 0 and dissipated when 
f{m) > 0. On a closed trajectory, the friction / must 
thus change sign. Since / is a decreasing function of to, 
this imposes /(O) > 0 for trajectories with m{z) > 0, or 
equivalently 


c > 




(C3) 


The conditions (Cl I 


and 


thus provide loose 


bounds on the subspace of the {c,Pg) plane which con¬ 
tains the three types of trajectories m{z) described above. 
These trajectories correspond to the three types of inho¬ 
mogeneous profiles seen in the microscopic models m- 
Before studying the stability and coarsening of these 
propagative solutions, we first need to understand pre¬ 
cisely how they are organised in the (c, pg) plane. In 
order to do so, we first analyze the phase portrait of the 
dynamics ( [T4| . We then study how the phase portrait 
evolves when pg and c are varied. 


As explained before, because of the condition (C2|, the 
extrema of at to = 0 , to_ , to+ are such that 0 < to_ < 
TO_|_, so that 0 and to_|_ are two maxima and to_ is a 
minimum of H. 

Linearizing the dynamics around one of the fixed 
points, we define to = toq -I- dm with toq = 0, to_,to+, 
so that TO = 8m and 


dz (<jm) " ^-H"{m^)/D -/(mo)/7?) (jto) 


The stability of the fixed points is given by the eigenval¬ 
ues Ai ^2 of the 2 x 2 matrix which read 


Ai,2(too) 


-fjrno) ^ // /(TOo) y 
2D \ 2D ) 


H"{mo) 

D 


( 22 ) 


• At the maxima toq = 0 and toq = m+ of El, 
H"(mQ) is negative and the two eigenvalues are 
thus real with opposite signs. These fixed points are 
saddle points with one unstable direction (Ai > 0) 
and one stable direction (A 2 < 0). 

• At the minimum toq = to_ of H, H"{mo) is posi¬ 
tive and the real part of the two eigenvalues have 
the same sign, given by —/(to_). The fixed point 
is stable when /(to_) > 0 and unstable when 
/(to_) < 0. Physically, when the friction of the 
Active particle is negative around to = to_, small 
perturbations are amplified, driving the trajectory 
away from the fixed point. Conversely, a positive 
friction damps any initial perturbation, leading to 
trajectories converging towards to_. 

When c and pg are such that f{m-) — 0, Ai _2 are 
complex conjugate imaginary numbers. A Hopf bi¬ 
furcation takes place, leading to the apparition of 
a limit cycle. 

At the onset of a Hopf bifurcation, a limit cycle ap¬ 
pears around the fixed point whose stability changes. In 
the following sections we elucidate how the interplay be¬ 
tween the saddle-point and the Hopf dynamics rules the 
non-linear dynamics of the Active particle and hence the 
polar-band shape. 


C. Stability of the fixed points 


The structure of the phase portrait is most easily cap¬ 
tured by locating the Axed points of (141 and studying 
their stability. We Arst rewrite (141 as a system of two 
Arst-order diAerential equations: 


d / m\ 
dz [mj 


/("») 

D 


H' (m) 


D 


( 20 ) 


D. Hopf bifurcation 

Let us Arst provide a comprehensive characterization 
of the Hopf bifurcation. It happens when the real part 
of Ai _2 vanishes, i.e., when 

f{m-) = c-^ - ^m-{c,pg) = 0 (23) 

where to_, which depends on both c and pg, is given by 
Eq. ([T7|. Equation (|23|) is satisAed on the line 


The Axed points are the solutions satisfying to = 0 and 
H'(m) = 0, i.e., the constant solutions extremizing H. 


Pgic) = Eg+ — 


-I- uoA)(—040^ -I- 04^0 A -I- vo^) 


(24) 















• When c > c*, the Hopf bifurcation is subcritical 
(A < 0). The system branches from an unstable 
fixed point when pg < (case IV, Fig. to an 
unstable limit cycle surroundiM a stable fixed point 
when Pg > p^ (case III, Fig. 

The organisation of these four typical cases in the (c, pg) 
plane is illustrated in Fig. We thus see that, when 
c < c*, limit cycles exist for pg smaller than p^, whereas 
when c > c*, they exist for pg larger than p^. The Hopf 
bifurcation line is thus a boundary of the domain of ex¬ 
istence of periodic propagative solutions of the hydro- 
dynamic equations. Let us now consider what happens 
when we explore the c, pg plane further away from the 
Hopf bifurcation line. 


E. Structure of the (c, pg) solution space 


FIG. 6. The four types of phase portrait (m, m) obtained 
in our system when changing pg and c. The line pg = pf 
is where the Hopf bifurcation takes place. The bifurcation is 
supercritical for c < c* and subcritical for c > c*. The plain 
(resp. open) black squares denote stable (resp. unstable) 
fixed points at m = m_. The plain (resp. dashed) black 
lines denote stable (resp. unstable) limit cycles. The crosses 
denote the saddle points at m = 0 and m = m+. The initial 
condition of each trajectory is marked by a magenta disk and 
the direction of “time” 2 ; indicated by a magenta arrow. 


which we call the Hopf transition line. 

Following standard text books in bifurcation the¬ 
ory m, the type of Hopf bifurcation (super- or sub¬ 
critical) is given by the sign of 


where uj = ^H"{m-, p^)/D > 0 is the imaginary part 
of the eigenvalues at the bifurcation point. Moreover 


dm_ —1 

-4a4(v?s - Pg) 


(26) 


is always negative because of condition The sign 

of A is thus given by the sign of Pg), which 

changes at c = c* with 


c* = \/^o(^'^4A-l-g) 
y/3a4 

Two different scenarios occur depending on whether c is 
larger or smaller than c*. 

• When c < c*, the Hopf bifurcation is supercritical 
(A > 0). The system branches from a stable fixed 
point for pg > p^ (case I, Fig. ^ to a stable limit 
cycle surroundingan unstable fixed point for pg < 
p^ (case H, Fig.j^. 


So far, we have shown that three different types of 
trajectories m{z) (periodic, homoclinic and heteroclinic) 
can be found by varying the values of c, pg. The sub¬ 
space where these physical solutions can be found was 
first bounded by the conditions (ICll), (|C2|), and ([C3|. In 


the previous section, we further found that the Hopf tran¬ 
sition line p^ (c) given by Eq. (241 is the upper boundary 
for the admissible values of pg when c < c* and the lower 
boundary when c> c*. 

To explore the remaining (c, pg) space, we numerically 
integrated the dynamical system (20 1 using a Runge- 


Kutta scheme of order 4. Starting from different initial 
conditions, one easily finds the basins of attraction of 
the different solutions. To locate unstable fixed points 
and limit cycles, we integrated the dynamics backward 
in time since they are attractors when 2 —>■ — 00 . As 
c and Pg vary, so do the shapes and sizes of the limit 
cycles. To quantify these variations, we measured the 
“amplitude” of a cycle, defined as the difference between 
the two extrema of m{z) 


Am = max[TO( 2 )] — min[m(z)] 


(28) 


We systematically vary pg at fixed c, first focussing 
on the case c < c* where the Hopf bifurcation is su¬ 
percritical. Decreasing pg, a stable limit cycle of van¬ 
ishing amplitude appears at pg = p^ (Fig. ^ panel A). 
The amplitude of the cycle then increases as pg decreases 
(Fig.[7l panel B) until it hits the fixed point at m = 0 
where the limit cycle becomes an homoclinic trajectory 
(Fig.0 panel C). For even lower pg the particle escapes to 
m = — 00 . The variation of the cycle amplitude with pg 
shown in Fig. can be qualitatively explained. When 
Pg decreases, the distance m_ — mj increases, where 
m/ = |(c — Auo/c) is the value of m where the fric¬ 
tion changes sign, i.e., /(to/) = 0. More energy is thus 
injected in the system and, to dissipate this energy, the 
trajectory need to go closer to to = 0 . 

A symmetric behavior is observed when c > c*, for the 
subcritical Hopf bifurcation. Increasing pg, an unstable 
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FIG. 7. Line at fixed c = 1.12 in the (c, pg) space, for c < c*. 
Top: Three phase portraits along the line, from the vicin¬ 
ity of the Hopf bifurcation line (A) until the homoclinic is 
found (C) Bottom-left: Size of the limit cycle Am defined 
in Eq. ( |28| . The limit cycle disappears when Am is large 
enough that the orbit reaches m = 0, where the trajectory is 
homoclinic. Bottom-right: Average density po of the solu¬ 
tions. Parameters: L> = wo = A = ^ = a 4 = (pg = l. 



FIG. 8. Unstable homoclinic trajectory starting from m = 
m+ for c = 1.16. Parameters: D = VQ = \ = ^ = a 4 = (pg = 

1 . 

limit cycle of vanishing amplitude appears at pg = . 

The amplitude of the cycle then increases with pg until 
the trajectory hits the point m = m+ where we have 
an (unstable) homoclinic solution that starts from m = 
TO_|_ as shown in Fig. [Sj Physically, when increasing pg, 
nif — rri- increases so^at the friction around the stable 
fixed point at m = m_ becomes larger and thus its basin 
of attraction (whose boundary is the unstable limit cycle, 
see Fig. becomes larger. 

All in all, the central results of this section is that 
all the admissible solutions lie in a band delimited by 
the Hopf bifurcation line p^(c) and a line where the ho¬ 
moclinic trajectories are found, as shown in Fig. In¬ 
side this band there exists stable non-degenerate limit cy¬ 
cles, corresponding to periodic propagating profiles. The 
unique heteroclinic trajectory is located exactly at c = c* 
where the Hopf bifurcation changes from supercritical to 
subcritical. We thus observe a 2-parameter family of pe¬ 
riodic solutions, a line of homoclinic trajectories and a 



FIG. 9. Space of all solutions. A 2-parameter family of pe¬ 
riodic orbits is found inside the band delimited by the Hopf 
bifurcation and the homoclinic trajectories. The constraints 
(G2) and (G3) are indicated by the black dashed lines. The 
constraint (Gl) lies out the range of this plot. Bottom: 
Zoom of the plot above around the point c = c* where 
the Hopf bifurcation changes from supercritical to subcriti¬ 
cal. This is also where the unique heteroclinic trajectory is 
found. The roman numbers refer to Fig. indicating the 
type of phase portrait found in each region. Parameters: 
D = Vo = A = 5 = 04 = yjg = 1. 

unique heteroclinic trajectory. Going back to the origi¬ 
nal pattern formation problem, they correspond to a 2- 
parameter family of micro-phase separated profiles, a line 
of isolated solitonic bands and a unique phase-separated 
state where a macroscopic polar liquid domain cruises 
through a disordered gas. 


F. Working at fixed average density 

In the microscopic models and the original hydrody¬ 
namic equations the average density po is a conserved 
quantity fixed by the initial condition. On the contrary, 
when considering the trajectories of the fictive particle 
m(z), po is not a priori fixed and varies between the 
different solutions. To compute the mean density on a 
trajectory m(z) we simply average p{z) = Pg + vom{z)/c 
over time. 

As shown in Fig. [^(bottom-right), we find that at fixed 
c < c*, Po decreases when pg decreases. It ranges from 
Po = Pg + when pg = p^ to po = Pg at the homo- 





































10 


clinic trajectory where the portion of the trajectory with 
m{z) yf 0 becomes negligibly small. Note that, at the het¬ 
eroclinic trajectory, pq can take a large range of values. 
Since the size of the gas and liquid domains are arbitrary, 
the average density can take any value in where 

Pg and p^ are the densities in the gas and liquid domains 
respectively. 

Fixing po adds a constraint that selects a line of solu¬ 
tions in the (c, Pg) space, as shown in Fig. [To] (left). For 
all po G these lines end at the heteroclinic tra¬ 

jectory. We also observe that, at fixed po, the closer the 
trajectories are to the heteroclinic solution, the larger 
their amplitude (see Fig. right). This means that 
along a line po = cst, the higher the amplitude the faster 
the band excitations propagate. This point will turn out 
to be crucial when discussing the coarsening dynamics at 
the hydrodynamic level in Sec. Em 

Until now, we have shown that three different types 
of possible trajectories m{z) exist, which correspond to 
all the propagative solutions observed in the microscopic 
models of flying spins. We have further identified the 
subset of values of the propagation speed c and the gas 
density pg for which these solutions exist. We can now 
turn to the study of their dynamical stability at the hy¬ 
drodynamic equation level. However, we first discuss an¬ 
alytically in the next section the shape of inhomogeneous 
solutions. 


G. Exact solution for the heterocline 



1.12 1.14 1.16 


FIG. 10. Top: Lines of solutions having a fixed average den¬ 
sity Po in the space of all solutions. Bottom: the amplitude 
Am of the cycles, defined in Eq. (281, along the lines po = cst 
increases when c increases. Parameters: D = vq = \ = ^ = 

a4 = (Pg = 1 . 


There are no general analytic solutions for the propa¬ 
gating inhomogeneous profiles. However, progress is pos¬ 
sible for some limiting cases. In the following we show 
that a complete solution for the heterocline - its position 
in the (c, pg) plane and its shape - can be determined 
exactly m- In 1^ we then show that, although exact 
solutions are not available, progress regarding the shape 
of the homoclinic solutions is achievable in the small D 
limit. 

To compute the shape of the heterocline, let us start 
from the ansatz 


7711,2(2) 


me 

~Y 


1 -I- tanh 


z — 



(29) 


Each of mi{z) and 7712 ( 2 ) describe an interface centered 
around 2 ; = 2 : 1,2 between a disordered phase with m = 
0 and an ordered phase with m = me. The complete 
heteroclinic trajectory then consists of two fronts glued 
together: An ascending front 777 . 1 ( 2 ;) with fci > 0 and a 
descending front 7712 ( 2 ) with ^2 < 0, with Z 2 S> 21 (see 
Fig. [TT| ); Being part of the same profile, the two fronts 
share the same celerity c and density pg. 

Moreover we know that me must be located at the 
second maximum of H so that 


^0 

2 a4C 




404 (yjg 

V 


2 

0 



(30) 


Plugging 7771 and m 2 in Eq. (131 and replacing me by its 


expression, one obtains for each of the fronts 


g{c, Pg, fci, 2 ) -I- h{c, Pg, fci, 2 ) tanh (^ 1 , 2(2 - 2 : 1 , 2 )) = 0 

(31) 

where g and h are complicated functions that we omit for 
conciseness. Eq. (31) can be true only if g and h vanish 


independently, for both ki and ^ 2 . 

We can express ki and k 2 as functions of c and pg 
by linearizing the ansatz (29) around m = 0. When 
^ 1 , 2(2 — 21 , 2 ) —t — 00 , one has 7711,2 ~ exp( 2 fci, 2(2 — 21 , 2 )) 
so that we can identify fci .2 with the two eigenvalues of 
the linear stability analysis Eq. (221. The ascending front 


is associated with the unstable direction ki = Ai/2 and 
the descending front with the stable direction ^2 = A 2 / 2 . 


Replacing fci ,2 by their values in Eq. (311, we have four 


equations for the two unknowns c and pg. After some 
algebra, one obtains a unique solution (c^, p^) with 


c = c = 


Pg=^9 


a/7;o(3o 4A -k ^) 

■s/SoI 

2i7o 


9714 A -(- 3^ 


(32) 

(33) 


This gives us the magnetization me and the the fronts 


me = m+ = 






















11 


0.6 - 



/ 


0.4 - 

/ 

1 


mi{z)l 

i 7712 ( 2 ;) 

0.2 - 

1 

1 


1 

1 

n n 

J 

V 

U.U 

z 

( 

1 

) 50 

1 1 

100 150 


FIG. 11. Comparison between the exact solution for the het- 
erocline (dashed lines) and the result from numerical integra¬ 
tion of Eq. (20 1 (blue line). Parameters: D = vo = X = ^ = 

a 4 = ipg = 1 . 


steepness ki ,2 as 


_ _ 2 uo__ 

\/ 3 a 4 t>o( 3 a 4 A -I- ^) 

^ _ y/vo{8a4D + C^) - 

4Z1 3CL 4 . (^30/4X ^) 

^ _ -\/vo{8a4D + g^) - yuOC 
4Zl'^3(i4(3(J4A ^) 


(34) 

(35) 

(36) 


In Fig. El we show that this solution matches exactly 
the heteroclinic orbit found by numerical integration of 
Eq. (pOl). 


V. BACK TO THE VICSEK-LIKE AND THE 
ACTIVE ISING MODELS: NON-LINEAR 
SOLUTIONS OF THE HYDRODYNAMIC 
EQUATIONS 


In Sec. |III| and |IV| we consider phenomenological hy¬ 
drodynamic equations and assumed the simplest possi¬ 
ble dependences of their coefficients with density. Here 
we extend our study to the more realistic hydrodynamic 
equations presented in Sec |III A| We first consider in 
Sec. |V A| the Vicsek hydrodynamic equations before turn¬ 
ing to the Ising hydrodynamic equations in Sec. |VB[ For 
sake of completeness, we also consider in Appendix [B] a 
more general case where the potential H appearing in 
the dynamical system for m{z) is not a polynomial in 
m. While not directly relevant for the hydrodynamic 
equations studied in this paper for the Vicsek and ac¬ 
tive Ising model, such Hamiltonians cannot be ruled out 
and may arise, for instance, from a density-dependence of 
the coefficient 04 in the phenomenological hydrodynamic 
equations 


A. Vicsek like equations 


Let us consider the hydrodynamic equations (3]|4) in¬ 
troduced by Bertin and co-workers to describe a simpli¬ 
fied version of the Vicsek model ng. These equations 
have the same structure as the phenomenological equa¬ 
tions (7 81 with two additional gradient terms |V(|toP) 
and —k(V • 'm)m. Importantly, all the coefficients 7 , v, 
K, and C. depend on the density. 

We follow the same approach as before. Looking for 
propagative solutions invariant in the transverse direc¬ 
tion we set rUy = 0 , write = m, and go to the comov¬ 


ing frame z = x — ct, to obtain 
Cp — VqTTI = 0 


(37) 


vfh + (c-— ymrh- -p + pTn — (^wX' = 0 

\ 2 c/ 2 

(38) 

Note that after setting niy = 0 the two k gradient terms 
of Eq. (|^ cancel each other. 

As before, Eq. (|3^ directly yields 


(39) 


P{Z) = Pg+ ^to(z) 


As we show in Appendix [C| Eq. (38 1 can also be greatly 
simplified using the explicit density-dependence of its co¬ 
efficients. Introducing 7 = 7 / 1 ^, C = C/^^i writing 
p = p-ip— p ,2 and = V 1 P+V 21 one obtains the second- 
order ordinary differential equation 

m -I- (a — ^m)m — a2m + asrri^ — 04771^ = 0 ( 40 ) 

where the coefficients are all function of c and Pg 


a = ( c - ^ ) {icipg + V 2 ) 


+ 7 ( 41 ) 


02 = PL2V2 + {P2V1 - k‘lV2)Pg - P-II'IP'I 

03 = — {2pgp.ll/i + P 1 V 2 - P2^l) 

C 


04 = c - 




Interestingly, Eq. (401 has exactly the same form as 
Eq. ( |13[ ) —the dynamical system stemming from the 
phenomenological hydrodynamic equations studied in 
Sec |IV] —although with coefficients depending in a more 
complicated way on pg and c. Their propagative solu¬ 
tions will thus be the same, up to a slightly different 
organisation in the c, pg plane. 

However, an important difference between the Vic¬ 
sek and phenomenological hydrodynamic equations, is 
the scaling of the magnetization with density in the or¬ 
dered phase. The Vicsek hydrodynamic equations (|3]|^ 
of Bertin and co-workers indeed predict that the homo¬ 
geneous ordered solution satisfies 


\mo\ 

Po 


1 

Po 


Po- 


PiVi 

c 


(42) 
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FIG. 12. Top left: Phase diagram for the propagative solu¬ 
tions of the hydrodynamic equations of Bertin et al. Eq. (40 1 . 
The color code and different phases are the same as for the 
04 = cst case Fig. We obtain the same three types of in¬ 
homogeneous profiles: periodic, homoclinic and heteroclinic, 
cr = 0.7. 


which is consistent with what is observed in microscopic 
models such as the Vicsek model. On the contrary, be¬ 
cause the coefficient 04 in Eq. ^ does not depend on den¬ 
sity, the phenomenological equations would yield 

^ 0 as po increases. In this region of parameter 

space, which is not the main focus of this paper, the Vic¬ 
sek hydrodynamic equations are thus more faithful to the 
phenomenology of microscopic models studied in Sec. [TT| 
than the phenomenological equations (7]|8l. As we show 
in Appendix however, these phenomenological equa¬ 
tions can recover a scaling akin to that of the Vicsek 
model if the coefficient 04 of Eq. ([^ is allowed to depend 
on density. This leads to a slightly more complicated dy¬ 
namical system m that we study in the appendix for 
completeness. 

Coming back to the dynamical system ( |40|, following 
the same method as that introduced in Sec. 13 one can 
derive analytical expressions for the Hopf bifurcation line 
p^{c) and the speed c* where the bifurcation becomes 
subcritical. We show in Fig.[^the phase diagram in the 
(c, pg)-plane and examples for the three types of inho¬ 
mogeneous trajectories. Again an exact solution for the 
heteroclinic trajectory can be derived and (the dashed 
lines in Fig. 121 is found at speed = c*. As expected, 
there is no qualitative difference with the simpler case 
studied in Sec. imi 


B. Active Ising model equations 

The active-Ising hydrodynamic equations m seem a 
priori different since they contain none of the non-linear 


gradient terms found in the phenomenological and Vicsek 
hydrodynamic equations. However, these terms are in 
fact generated by the diffusion term in the dynamical 
equation (§ for the density field m, so that we will once 
again recover the same types of propagative solutions. 

Looking for stationary profiles in the comoving frame 
z = X — ct, Eq. © and (§ reduce to 

Dp + cp — vom = 0 (43) 

/ 7* \ 

Dm + cm — vqp -1-2 /3 — 1-|m — = 0 (44) 

V PJ P 


Equation (43) can be solved by expanding p in gradients 
of TO. Introducing the ansatz 


p{z) = pg + Y^ 


Oik 


k=0 


d^m 

dz^ 


(45) 


into Eq. (|43|), and solving order by order, we get the 

D 


following recursion relation 
Vo 


0-0 — 


from which we obtain 


Clk+l — 




k=0 


D 


-ak 




dz^ 


(46) 

(47) 


In the following we retain only the first order terms in 
gradient: 


, , Vo , . Dvo ... D^vo .. . . 

P\D = Pg + —rn[z) - ^rn[z) + —^m{z) 


(48) 


To simplify Eq. (441 we linearize around the density 
ipg = r/(/3 — 1) where the disordered profile becomes 
linearly unstable. As shown in m, this is a good ap¬ 
proximation when /3 > 1 because ipg ^ 00 while p — pg 
remains finite for inhomogeneous solutions. We then ob¬ 
tain 

2 /’ TO^ 

Dm + cm — Vop-\ - k(p — Pn)m — a^r = 0 (49) 

^9 ‘Pg 


Finally, inserting Eq. (48 1 in the previous equation 
we obtain a second-order ordinary differential equation 
with exactly the same terms as those obtained from the 
phenomenological and Vicsek hydrodynamic equations, 
Eq. (pi) and (|40|: 




Drh + 
with 

D = D 




^m\ m — a^m + azm?— apm^ = Q (50) 


1 + 4 




02 = 


2 r((/ 3 g - Pg) 




03 = 


ArDvo 

(c2 -h vDpI 


2rvo 

cpl 


04 = 


(51) 

(52) 


This equation has again the same qualitative behavior 
as in the phenomenological theory and the same three 
types of inhomogeneous solutions are found with the 
same organization in the (c, pg) parameter space. 
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VI. LINEAR STABILITY OF THE 
PROPAGATIVE SOLUTIONS IN THE ID 
HYDRODYNAMIC EQUATIONS 

In sections |IV| and |V] we have found and classified 
all the propagative solutions of different hydrodynamic 
equations. We have shown that in all cases three types 
of such solutions exist periodic patterns of finite-size 
bands, solitary band solutions, and phase-separated so¬ 
lutions. These solutions were found as stable limit cycles, 
homoclinic, and heteroclinic orbits m{z) of the reduced 
dynamical system (13l. This study does not tell us any¬ 


thing about the local and a fortiori global stability of 
these solutions as solutions m{x — ct) and p{x — ct) of 
the original hydrodynamic partial differential equations. 
Indeed, a stable orbit of the reduced dynamical system 
can very well be unstable to spatiotemporal perturba¬ 
tions when re-expressed as an inhomogeneous propaga¬ 
tive solution of the hydrodynamic equations. For ex¬ 
ample, consider the fixed point at m = m_ which is 
stable for pg > in the region I and III in Fig. ^ 
Without performing any calculation, we know that the 
corresponding homogeneous solution (po = Pg + 
mo = m_) is unstable in the hydrodynamic equations 
because it lies inside the spinodal region where no homo¬ 
geneous solution is linearly stable. 

In this section we study the linear (local) stability of 
these solutions at the hydrodynamic level. Although, it is 
a well-defined linear problem, determining the stability of 
inhomogeneous solutions of nonlinear partial differential 
equations cannot be done analytically even in the (rare) 
cases where these solutions are known analytically. A di¬ 
rect numerical study is possible (for an example, see, e.g. 
|25|'). but is rather tedious, all the more so as we deal with 
2D hydrodynamic equations. In the following, we study 
mostly the linear stability of the hydrodynamic equations 
reduced to one dimension (that of propagation), using a 
simple numerical procedure explained below. For an ac¬ 
count of some fully 2D preliminary investigation, see the 
end of this section. 



FIG. 13. Stability of homoclinic orbits in the hydrodynamic 
equations. Auc is the difference in the amplitude of the 
solution measured between t = 0 and t = Tg. The solu¬ 
tion is declared to be linearly stable if Sm{Ts), defined in 
Eq. (531 is smaller than a threshold that we take at 10“®. 


uo = A = 5 = 04 = Pg = 1, c = 1.125 (right) 


and periodic boundary conditions. 

Preparing the initial condition of the hydrodynamic 
equations always brings in discretization errors. Indeed, 
because of the periodic boundary conditions in the hy¬ 
drodynamic equations, we need to select a portion of the 
solution m{z) which is a multiple of the period. This is 
done with an error of order dx, the space discretization 
step used in the numerical integration of ([^ and (|^. 

Accordingly, we observe a rapid relaxation at short 
times due to the discretization errors. Subsequently, we 
find that the original solution is either stable at long 
times or is quickly destabilized and converges to another 
solution. To analyze systematically the stability of the 
propagative solutions, we defined a quantitative criterion 
for the stability: We choose a time Tg = 2000 much 
larger than the relaxation time of the initial perturba¬ 
tion (which happens in a time ~ 100) but not too large 
to test solely the linear stability regardless of a possible 
long-time coarsening dynamics that could be induced by 
numerical noise. We then measure the amplitude of the 
solution I Am | (t) defined in Eq. (28) as a function of time. 
If 


5m{Ts) = \Am{Tg) - Am(0)| 


(53) 


A. Numerical procedure 

We investigated numerically the stability of the prop¬ 
agative solutions of the phenomenological hydrodynamic 
equations Eq. (7][8l reduced to one space dimension, that 
of propagation. To do so, we select a solution of the corre¬ 
sponding classical mechanics problem. (13) and use it as 
initial condition of the numerical integration of Eqs. 0 
and 

The numerical integration is done using a semi-spectral 
algorithm, the linear terms being computed in Fourier 
space, the non-linear ones in real space and a semi- 
implicit time-stepping. This method, where the fields 
p and m are represented by their N first Fourier modes 
(with N large enough that the simulation has converged), 
is well suited to simulate systems with diffusion terms 


is smaller than 10 the solution is said to be stable 


and unstable otherwise (see Fig. 131. This protocol does 
not give exact answers to the question of linear stability, 
since in particular the small but finite initial perturba¬ 
tions may take the initial condition out of the basin of at¬ 
traction of a (stable) solution. But the results presented 
below are relatively robust to changing our numerical res¬ 
olution and the conditions used to decide stability, and 
we are thus confident that they represent well the ’true’ 
subset of linear stable solutions. 

The precise criterion does not matter much for the re¬ 
sults since we find an abrupt transition from stable to 
unstable trajectories (see Fig. 131 which is visible on a 
variety of observables (norm of the fields, period, max 
and min values, etc). 
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B. Results 

Figure [^contains one of the central results of our pa¬ 
per: only a very small subset of the propagative solutions 
are stable at the hydrodynamic level. However this sub¬ 
set includes the three possible types of trajectories: peri¬ 
odic bands, solitonic bands, and phase-separated profiles. 
The linearly stable solutions are all found in the region of 
the (c, Pg) plane close to the heteroclinic trajectory and 
close to the line of homoclinic orbits. Examples of stable 
and unstable solutions are shown in Fig. 

To understand why only the vicinity of the heteroclinic 
and homoclinic solutions is stable, one can argue that the 
solutions must have a large enough amplitude to be dy¬ 
namically stable. Indeed, the periodic solutions oscillate 
around m = m_ which lies inside the spinodal region 
of the hydrodynamic equations, where no homogeneous 
solutions is stable. Small-amplitude oscillations around 
this point should thus also by dynamically unstable and 
only large enough amplitude cycles, found near the ho¬ 
moclinic line and the heteroclinic trajectories, are stable. 

Note that the region where stable solutions are found 
has a rather rough boundary in Fig. This is most 
probably an artefact due to the initial discretization er¬ 
ror which is not controlled and varies from one propaga¬ 
tive solution to another. Close to the threshold of linear 
instability, this can easily make a (linearly) stable trajec¬ 
tory unstable. 

We report finally on fully 2D simulations performed on 
rectangular boxes of width Ly — 100, with small noise 
added to each grid point on a ID solution extended triv¬ 
ially along y. These yielded essentially no change with 
respect to the results presented above. While this encour¬ 
ages us to believe that no unstable mode has components 
along y, and thus that the subset determined above cor¬ 
responds to the linearly stable solutions of the full 2D 
hydrodynamic equations, we remain cautious. As a mat¬ 
ter of fact, recent results obtained in the case of hydro- 
dynamic equations for active nematics have revealed the 
existence of (very) long wavelength instability along the 
homogeneous direction of band solutions. m A similar 
investigation is left for future studies. 


VII. COARSENING IN THE HYDRODYNAMIC 
EQUATIONS 

The two-dimensional subset of propagative solutions 
that are linearly stable still contains an infinity of so¬ 
lutions, including smectic patterns, solitary bands, and 
phase-separated profiles. We can thus wonder whether a 
single solution is selected in the hydrodynamic equation 
starting from a random initial condition. 

One can obtain some insight into this question by 
studying the lines of propagative solutions having a given 
average density poi shown in Fig. [T^ Indeed, po is con¬ 
served in the hydrodynamic equations and any stable 
propagative solutions has to lie on such a line. As dis¬ 



FIG. 14. The blue region is the subset of admissible prop¬ 
agative solutions which are linearly stable. None of the so¬ 
lutions for c > c* is stable in the hydrodynamic equations. 
Do = A = ^ = 04 = Pg = 1. We use system sizes = 300 
or more (adapted to fit the solution) with resolution da; = 0.5 
and dt = 0.1. 



FIG. 15. Evolution of two propagative solutions in the hy¬ 
drodynamic equations The initial conditions are sta¬ 

ble limit cycles of the dynamical system (13l for m(z). Top 
row: unstable solution with c = 1.135, pg = 0.8351. Bottom 
row: stable solutions with c = 1.14, pg = 0.8341. System size 
340 X 100, dt = 0.1, dx = 0.4. uq = A = ^ = a4 = pg = 1. 


cussed earlier, the larger the amplitude of a solution, the 
faster the propagation. We can thus expect that if several 
travelling bands coexist, they will encounter and merge 
because of their different speeds. This will in turn in¬ 
crease the sizes of the surviving objects and hence their 
speed. This mechanism would naturally lead the sys¬ 
tem toward the heteroclinic solution, i.e., to the phase- 
separated state. 

We checked this scenario by simulating the phe¬ 
nomenological hydrodynamic equations (7][8) with three 


different initial conditions, as shown in Fig. 16 


• We build an initial condition made of two homo- 
clines glued together, both solutions of the differ¬ 
ential equation (13) for different speed c and inside 


the stability domain of Fig. To avoid discon- 
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tinuities we interpolate smoothly between the gas 
densities of the two solutions using a hyperbolic 
tangent function. We then observe that the two 
travelling bands get closer until, when close enough 
(though not in contact), the smaller one evaporates, 
its mass being transfered to the second band. The 
final solution is hence a single larger isolated band. 


• Starting from a random ordered solution with con¬ 
stant p and a magnetization m fluctuating around 
mo yf 0, several bands form at short times that all 
go in the same direction. The system then coarsens 
because of the speed differences between the liq¬ 
uid droplets until only “phase-separated domains”, 
that all have the same speed, remain. In practice, 
because the band speeds can be very close, a fi¬ 
nal state with a single phase-separated profile may 
not be reached within the time-scales of our simula¬ 
tions and a precise study of this coarsening regime 
is beyond our numerical capacities. 


• Starting from a disordered initial condition with 
in{x) fluctuating around 0 and a density inside 
the spinodal region, many domains of positive and 
negative magnetization form. These objects then 
encounter and merge yielding a rapid coalescence 
process. Because of the periodic boundary condi¬ 
tions, this process typically yields a single phase- 
separated state for the system sizes we considered. 
For larger systems, the coalescence process could 
result in a larger number of bands propagating in 
the same direction. We should then observe the 
same type of coarsening as in the second case dis¬ 
cussed above. 


At the level of hydrodynamic equations, the fact that 
larger ordered domains travel faster leads to a natu¬ 
ral coarsening towards the phase-separated states. This 
coarsening relies both on a coalescence process and, when 
bands travelling in the same direction are close enough, 
on a ripening during which a smaller band evaporates 
and is “swallowed” by its larger neighbour. 


Note that we studied the stability of propagative solu¬ 
tions and their coarsening process using the phenomeno¬ 
logical hydrodynamic equations (7][8). While the precise 
results will depend on the set of hydrodynamic equations 
under study, our simulations of the Vicsek and Ising hy¬ 
drodynamic equations, (3]|4| and (^), did not suggest 
any qualitative difference. Their comprehensive investi¬ 
gation, which is beyond the scope of this paper, would 
nevertheless be interesting, especially to quantifiy the 
role played by the non-linearities in the vectorial hydro- 
dynamic equations Q. 


VIII. PHASE DIAGRAM OF 
COARSE-GRAINED HYDRODYNAMIC 
EQUATIONS 


We can now construct the full phase diagram for 
the Ising and Vicsek hydrodynamic equations in the 
temperature-density ensemble (or noise-density for Vic¬ 
sek), Fig. 17 It is composed of two spinodal lines, which 
are the limit of linear stability of the homogeneous disor¬ 
dered and ordered profiles, and two binodal lines outside 
which inhomogeneous solutions cannot be observed. 


A. Spinodal lines 


The spinodal lines can be determined analytically from 
a linear stability analysis, as was previously done for the 
Vicsek model [IS] and the active Ising model [13]. The 
lower spinodal line ipg is simply the density at which the 
coefficient of the term linear in m in the hydrodynamic 
equations changes sign (in our phenomenological equa¬ 
tion, when 02 = 0). It reads for the Ising and Vicsek 
hydrodynamic equations 




V ( \ A*2 


7r(l - e-'^'/^) 

Pi ~ 4(e-‘"V2 _ 2/3) 


(54) 

(55) 


where [3 = 1/T. The higher spinodal Lpi can be deter¬ 
mined exactly for the Ising hydrodynamic equations 


T’l {T) = (fig 


Vo 


a {vqK + 8D(l3 — 1)^) -I- VqK + 8Da{(3 — I) 


2vqK + 8Da{j3 — I) 


(56) 

where k = 2 + a — 2f3 and a = /3^(1 — /3/3). For the 
Vicsek hydrodynamic equation, the exact determination 
of (p£ is much more cumbersome m- In Fig.nUwe show 
the line ipY (a) computed numerically by simulating the 
Vicsek hydrodynamic equations at different densities in 
the homogeneous ordered state and looking at the growth 
of a small perturbation. 


B. Binodal lines 

The binodal lines pg and pe are defined as the minimum 
and maximum global densities beyond which inhomoge¬ 
neous propagative profiles cannot be observed in sim¬ 
ulations of the hydrodynamic equations. As explained 
in Sec. HVFj at the heteroclinic trajectory the size of 
the liquid and gas domains are arbitrary so that phase- 
separated solutions can have any density in the range 
[Pg^Pe]- that, for all parameters we tested. 

Pi = Pe, i-e., no other stable solution has a larger density 
than the liquid domain of the heteroclinic solution. 

The situation is more subtle for the lower binodal pg. 
Depending on the external parameters, the line of ho¬ 
moclinic trajectories solution of the dynamical system 
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FIG. 16. Coarsening in the phenomenological hydrodynamic equations (^ 8 l from three different initial conditions. Top: Two 
homoclinic trajectories glued together. Middle: Disordered initial condition. Bottom: Ordered initial condition. System 
size 500 X 100 (top and bottom), 1000 x 100 (middle). The profiles are averaged along the y-direction in which the system is 
invariant, wq = A = ^ = 04 = ipg = 1 , dt = 0 . 1 , dx = 0 . 1 . 



FIG. 17. Phase diagram of the Visek (left) and Ising (right) hydrodynamic equations. The spinodal lines yjg and (fie are known 
exactly except 93 ^ in the Vicsek hydrodynamic equation which we computed numerically by looking at the stability of systems 
of size 50 X 50. The binodals pg and pi are the coexisting densities of phase-separated (heteroclinic) solutions as explained 
in the text. The dashed line indicate the asymptote above which only disordered solutions exist. Insets are close-ups on the 
high-density regions with a logarithmic scale on the a;-axis. Parameters: i!o = D= l = r = l (Ising), vo = 1 (Vicsek). 
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0.82 0.84 0.86 0.88 0.76 0.78 0.80 


FIG. 18. Propagative solutions of the Vicsek hydrodynamic 
equations for a = 0.2 and 0.7. The solution with the lower 
mean density pmin is found at the minimum of the line of 
homoclinic solutions. At low noise, it does not coincide with 
the heteroclinic solution (blue squares), uq = 1 . 


is not always monotonous as a function of c. For ex¬ 
ample, in the Vicsek hydrodynamic equation it is not 
monotonous at low noise, as shown in Fig. [T^ In this 
case the minimum density accessible to propagative so¬ 
lutions is the minimum of the homoclinic line. This so¬ 
lution need not be stable in the hydrodynamic equations 
so that one should repeat the stability analysis done in 
Sec. I VI I for each value of the noise to determine exactly 
the binodal line. For simplicity, the lower binodal of the 
Vicsek hydrodynamic equation shown in Fig. is the 
gas density read from the phase-separated profile, which 
is true at high noise and a good approximation of pg at 
lower noise values. 

For the Vicsek hydrodynamic equations, the coexisting 
densities of the heteroclinic trajectory are known exactly 
whereas in the Ising case they can be determined analyt¬ 
ically only after linearizing Eq. (441 around p = (pg. The 
binodal lines in the phase diagram of the Ising hydrody¬ 
namic equations shown in Fig.j^are thus determined by 
integrating numerically the hydrodynamic equations and 
measuring the density of the liquid and gas domain of a 
phase-separated solution whereas we plot the analytical 
solution in the Vicsek case. 


IX. CONCLUSION 


Before summarizing our results, let us discuss how the 
study of the hydrodynamic equations presented in this 
article compares with the phenomenology of the micro¬ 
scopic models. The phase diagrams of the Ising and Vic¬ 


sek hydrodynamic equations shown in Fig. 17 are quali¬ 
tatively similar. They are also consistent with the phase 
diagrams of the microscopic models shown in Fig. The 
hydrodynamic equations thus provide a picture which 
is consistent with the liquid-gas framework discussed in 
Sec. lie For instance, the asymptote at finite noise (or 


temperature) can again be seen as the simplest way of 
forcing the system to cross a transition line to go from 
its gas to liquid phase. 

The comparison between the phase-separated regions 
of the microscopic models and hydrodynamic equations 


is however more subtle. We have indeed shown, using 
the phenomenological hydrodynamic equation 8), that 
the coarsening leads, at the hydrodynamic level, to the 
phase-separated solution. This is true in the active Ising 
model but not in the Vicsek model, which only shows 
micro-phase separation and thus a coarsening leading 
to a periodic solution (see Fig. [^. This difference be¬ 
tween microscopic models and hydrodynamic equations 
is however not surprising since it was recently shown that 
fluctuations are essential to account for the selection of 
micro-phase separated profiles in the Vicsek model [Hi- 


More precisely, phase-separated and micro-phase- 
separated solutions are both linearly stable in the hy¬ 
drodynamic equations for vectorial order parameters. As 
noise is added to these equations, though, large bands 
are destabilized and break in periodic array of finite-size 
bands, in agreement with what is seen in microscopic 
simuations of the Vicsek model [Tl]. We can now ten¬ 
tatively put together these results to propose a coarsen¬ 
ing process that would lead to the micro-phase separated 
states seen in the Vicsek model. Starting from a profile 
with many coexisting bands, the coarsening at fixed po 
would tend to lead to the formation of larger and larger 
bands. This coarsening would however be arrested by 
the fact that the fluctuations in the Vicsek model set a 
maximal size beyond which bands are (non-linearly) un¬ 
stable. How this size is selected however remains to be 
determined. 


All in all, we have shown in this paper that the hy¬ 
drodynamic equations describing polar flocking models 
generically share the same propagative solutions. We 
found three types of solutions that are all observed in 
microscopic models: periodic orbits that describe mi¬ 
crophase separation, homoclinic orbits describing soli- 
tonic objects and heteroclinic trajectories that corre¬ 
spond to phase separation. Only a subset is stable in the 
hydrodynamic equations, but it still contains the three 
types of solutions. 

The coarsening in the hydrodynamic equations favors 
the fastest solutions which are also the largest patterns. 
It thus leads naturally to the phase-separated solution 
which controls the phase diagram of the hydrodynamic 
equations. The same behavior is observed in the micro¬ 
scopic active Ising model whereas one can understand 
the microphase separation of the Vicsek model as an ar¬ 
rested coarsening, the largest pattern being non-linearly 
unstable to the (giant) fluctuations of the liquid phase. 
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Appendix A: Limit of small D 


While we do not have an analytic solution for the ho¬ 
moclinic profiles, lots of insight on their physics can be 
gained by studying the limit of small diffusion coefficient 
D. This is most easily done by introducing the auxiliary 
variable Z = Dm -I- F(rn) where 


F{m) 



(Al) 


such that F'(m) = /(m). Our dynamical system can 
then be recast as 


dz{z)-[ -H'{m) ) 


(A2) 


Let us consider a large-amplitude orbits that start close 
to m = 0. An example of such orbits and the correspond¬ 
ing phase portrait for the (m, Z) variables is shown in 
Fig. When D is small, m relaxes quickly to zero so 
that Z relaxes to the parabola Z = F{m) in a time ^ D. 
Following the trajectory in Fig. starting from point 
A at m « 0, the trajectory between A and B is above 
the parabola F{m). The distance with the parabola first 
increases when m < to_, which implies F['(m) < 0 and 
Z > 0, and decreases afterwards when m > m_. The 
distance with the parabola stays of order D, set by the 
relaxation time of m. When Z = /(m) at point B, m 
relaxes to m = 0 in a time ^ D. This gives profiles with 
sharp leading fronts and long exponential tails, which are 
indeed consistent with the profiles seen at small temper¬ 
ature in the Vicsek model and its putative hydrodynamic 
description nmsiiiH]. 

This picture is consistent with the fact that at leading 
order in D the eigenvalues at point (m = 0, m = 0) read 


Pg 

c — Xvq jc 


^2 


c — Auo jc 
D 


(A3) 


so that we have a slow unstable direction Ai and a fast 
stable direction A2. These two eigenvalues indeed control 
a large part of the trajectory, as shown in the right panel 
of Fig. [1^ 


Appendix B: Phenomenological equations with a 4 (p) 


In the hydrodynamic equation (|^ the homogeneous 
ordered solutions have a magnetization \fho\ = \/ 

Thus, in the simplified case studied in Sec. |TTl] where 
02 = p — and 04 = cst, one observes that 


l^ol ^ j p-^g _ 

p V p->-oo 


(Bl) 


On the contrary we observe in the microscopic Vicsek 
and active Ising models that, at large densities, |mo|/po 
reaches a constant Po < 1 (set by the noise intensity 



FIG. 19. Left: Phase portrait for the variables m and 
Z = Dm+F{m) where F{m) is the red parabola. Right: The 
corresponding trajectory. The dashed lines are the linear ap¬ 
proximation around m = 0. The inset in semi-log scale shows 
the exponential tail close to m = 0. Parameters: D = 0.01, 
no = A = 5 = 04 = P 9 = 1. 


in the microscopic models). This can be achieved by 
assuming that 04 = (pPp )“^ as in [T5] so that 


\rno\ 

P 



- 

p^oo 


(B2) 


Repeating the same treatment as before, we look for 
Id propagative solutions with speed c which are the so¬ 
lutions of 


Dm + {c -^ — ^m)fn — {(pg — Pg)m 

W 2 _ 

c™ Po{Pg + ^rn) 


= 0 (B3) 


which can be written in the same form as Eq. (141 with 
a different potential 


Dm = —f{m)m — F{' (m) 


(B4) 




C^Pg 


(y’g - Pg) 


2 P2 


0 


3P2u 


Wo 

3c 


cVg 

P^vl 


f{m) = c — 

c 


log{vom + cpg) (B5) 
(B6) 


We find a behavior very similar to the case 04 = cst 
discussed in section IV Under the same constraints (Cl- 
3), the potential H(m) has the same form with maxima 
in TO = 0 and to = to+ and a minimum at to = to_ given 

by 


m± = 


CWoPq (2pg - (fig) ± CPoA/(PoWoPg)2 - 4:C^Pg{(Pg ” Pg) 
2(c2 - P>2) 


(B7) 

One observes the same three types of trajectories al¬ 
ready shown in Fig. periodic, homoclinic and hetero- 
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IS 


clinic. The phase diagram (pg, c), shown in Fig. 
also similar to the case 04 = cst except in a very small re- 
gion close to the heteroclinic trajectory. We find again a 
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dynamical system as for the phenomenological equations, 
but with coefficients depending in a more complicated 
way on pg and c. 


We start from the hydrodynamic equations (3]|4| in the 
comoving frame z = x — ct: 


Cp — VqTTI = 0 


(Cl) 


2c 


^^0 . 


vfh + [c — — \ fn — 'ymm — 7 rP + P’m- — = 0 


(C2) 


The coefficients in Eq. (C21 are greatly simplified by di¬ 
viding the equation by u, thus getting 


2c 


m — ^rnrh——p+—rn — (rn^ = 0 (C3) 
2 v V 


FIG. 20. Top row: Phase diagram for the propagative so¬ 
lutions of Eq. (B4|. The right plot is a zoom on the region 
c* < c < c^. Phase I to IV are the same as in the case 04 = cst 
(see Fig. |^. Bottom row: Phase portraits of the solutions 
in phase V for the same c = 1.2421 and pg = 0.81228 (left) 
and pg = 0.81248 (right). When pg increases, the size of the 
stable cycle decreases while the size of the unstable cycle in¬ 
creases, until they collide as indicated on the phase diagram. 
Vo = A = ^ = Po = ipg = 1. 


where 7 = ^jv and C^ = C^jv depend only on the noise a 
on the alignment interaction m 
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(C4) 

(C5) 


2-parameter family of periodic orbits and a line of homo¬ 
clinic solutions which terminates at a unique heteroclinic 
trajectory. 

As before, we can compute analytically the line p^(c) 
where the Hopf bifurcation takes place and the speed c* 
such that the bifurcation is supercritical for c < c* and 
subcritical for c > c*. The difference with the previ¬ 
ous case is that the heteroclinic trajectory, for which we 
do not have an exact solution anymore, is not found at 
c = c* but at > c*. As a consequence, we observe a 
new phase (shown as number V in Fig. 201 in which two 
limit cycles are found. A large stable cycle surrounds a 
smaller unstable limit cycle which itself encapsulates the 
stable fixed point m = m_. As seen from the phase por¬ 
traits in Fig. 20 when increasing pg the size of the stable 
cycle decreases whereas the size of the unstable one in¬ 
creases. Thus the upper limit of existence of this phase 
(the magenta line in the phase diagram) is the moment 
where the cycles collide and annihilate. 


Appendix C: Propagative solutions of the Vicsek 
hydrodynamic equations 


while p, and iz depend also on the density 




^-1 = 


Stt 




+ (1 - e--^) 


(C6) 

(C7) 


One can solve Eq. ( |C1[) to get p{z) = Pg + '^m{z) as 
before and use it in Eg. ( [C3| ) to obtain the second-order 
ordinary differential equation 


m -I- (a — ^m)rn — a 2 m + = 0 (C8) 

where the coefficients are all function of c and pg 

a= {c-Y^{vipg + V2) 

02 = P 2 V 2 + {P2Vl - Pl^2)Pg - PiVipI 

03 = — { 2 pgPlVl + P1V2 - P2V1) 

c 

p '^0 

04 = C-ypio 


Here we show how the study of propagative solutions 
in the Vicsek hydrodynamic equations leads to the same 


where pi ,2 are defined hy p = pip — P 2 and 121^2 by iz ^ = 
l^lP+ V2- 


[1] Sepulveda, N., Petitjean, L., Cochet, O., Grasland- Biol., 9(3), el002944 (2013) 

Mongrain, E., Silberzan, P., Hakim, V. PLoS Gomp. 



















20 


[2] F. Peruani, J. Starruss, V. Jakovljevic, L. Sogaard- 
Andersen, A. Deutsch, M. Bar, Phys. Rev. Lett. 108, 
098102 (2012) 

[3] V. Schaller, C. Weber, C. Semmrich, E. Frey, A. R. 
Bausch, Nature 467, 73-77 (2010). 

[4] Y. Sumino, K. H. Nagai, Y. Shitaka, D. Tanaka, K. 
Yoshikawa, H. Chate, K. Oiwa, Nature 483 448-452 
( 2012 ). 

[5] J. Deseigne, O. Dauchot, H. Chate. Phys. Rev. Lett., 
105, 098001 (2010); J. Deseigne, et al. Soft Matter 8, 
5629 (2012); C.A. Weber, et al. Phys. Rev. Lett., 110, 
208001 (2013). 

[6] A. Bricard, J.B. Caussin, N. Desreumaux, O. Dauchot, 
and D. Bartolo, Nature, 503, 95 (2013) 

[7] Shashi Thutupalli et al 2011 New J. Phys. 13 073021 

[8] T. Sanchez; D. T. N. Chen, S. J. DeCamp, M. Heymann, 
Z. Dogic, Nature 491 431 (2012). 

[9] T. Vicsek, A. Czirok, E. Ben-Jacob, 1. Cohen, and O. 
Shochet, Phys. Rev. Lett. 75, 1226 (1995) 

[10] Mermin, N. D., Wagner, H. Phys. Rev. Lett. 17, 1133 
(1966) 

[11] J. Toner, and Y. Tu, Phys. Rev. Lett. 75, 4326 (1995). 
J. Toner, Y. Tu, and M. Ulm, Phys. Rev. Lett. 80, 4819 
(1998) 

[12] G. Gregoire, and H. Chate, Phys. Rev. Lett. 92, 025702 
(2004) ; H. Chate, F. Ginelli, G. Gregoire, and F. Ray¬ 
naud, Phys. Rev. E 77, 046113 (2008) 

[13] A.P. Solon, and J. Tailleur, Phys. Rev. Lett. Ill, 078101 
(2013) 

[14] A. P. Solon, H. Chate, J. Tailleur, Phys. Rev. Lett. 114, 
068101 (2015) 


[15] E. Berlin, M. Droz, G. Gregoire, Phys. Rev. E 74 022101 
(2006); J. Phys. A 42 445001 (2009); A. Peshkov, E. 
Berlin, F. Ginelli and H. Chate, Eur. Phys. J Special 
Topics 223, 1315 (2014). 

[16] J-B Caussin, A. Solon, A. Peshkov, H. Chate, T. Dauxois, 
J. Tailleur, V. Vitelli, and D. Bartolo, Phys. Rev. Lett. 
112, 148102 (2014) 

[17] S. Mishra, A. Baskaran, M.C. Marchetti. Phys. Rev. E 
81, 061916. (2010); A. Gopinath, M.F. Hagan, M.C. 
Marchetti, A. Baskaran. Phys. Rev. E 85, 061903. (2012) 

[18] T. Ihle, Phys. Rev. E 83, 030901 (2011); T. Ihle, Phys. 
Rev. E 88, 040303(R) (2013). 

[19] F. D. C. Farrell, J. Tailleur, D. Marenduzzo, M. C. 
Marchetti, Phys. Rev. Lett. 108, 248101 (2012) 

[20] R. Grofimann, L. Schimansky-Geier, and P. Romanczuk, 
New Journal of Physics 15, 085014 (2013) 

[21] Y. A. Kuznetsov, Elements of Applied Bifurcation The¬ 
ory, Second Edition. Springer. (2004) 

[22] In the Vicsek model, the region in which solitary bands 
can be seen vanishes as the system size is taken to infinity. 
At fixed intensive parameters, doubling the system size 
indeed doubles the number of bands and increasing the 
system size hence asymptotically leads, in principle, to a 
smectic arrangements of bands. 

[23] Note that m is actually a momemtum field in this con¬ 
text. 

[24] Again, in the infinite-size limit, the homoclinic or¬ 
bits/solitary wave solutions are only observable in a codi¬ 
mension one subset of parameter space. 

[25] H. Chate, and P. Manneville, Phys. Lett. A 171, 183 
(1992). 

[26] S. Ngo, et al. Phys. Rev. Lett. 113, 038302 (2014). 



